// gen_soft.cc
// Generate one soft low pT event 
// Used for merging with toyMC events
// 1 October 2011

#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <time.h>
#include <getopt.h>
#include <cstdlib>
#include <iomanip>

#include <TROOT.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TF2.h>
#include "tppmc/TF2J.h"
#include "tppmc/TF1J.h"
#include <TFile.h>
#include <TObject.h>
#include <TNtuple.h>
#include <TLorentzVector.h>
#include <TTree.h>
#include <TRandom.h>
#include <TRandom3.h>
#include <TClonesArray.h>
#include <TGenPhaseSpace.h>

#include <sstream>
#include <vector>

#include <TMath.h>
#include <TPythia6.h>
#include <TMCParticle.h>
#include <TVector3.h>

#include "tppmc/TAcceptParticle.hh"

using namespace std; 

#define RHIC200_CMS_ENERGY          200.0
#define RHIC500_CMS_ENERGY          500.0

#define pytune pytune_
extern "C" int  pytune(int *itune);

TRandom *rang = new TRandom3(static_cast<int>(time(NULL)));

typedef struct{
  double px; 
  double py;
  double pz;
  double e;  
  int type;
} pythpart; 



// MAIN PROGRAM
//--------------------------------

int main() {//int argc, char **argv) {

  TPythia6 EvtGen;

  // QCD low pT physics process
  EvtGen.SetMSEL(0);
  EvtGen.SetMSUB(95,true);

  // Pythia TuneA:
  int tuneval = 100; 
  pytune(&tuneval);

  // Set random seed
  time_t t1 = time(NULL); 
  EvtGen.SetMRPY(1, static_cast<int>(t1));     //random seed number

  EvtGen.Initialize("CMS","p","p",RHIC200_CMS_ENERGY);

  // Generate a pythia event, get particles

  EvtGen.GenerateEvent();
  TClonesArray *particleArray =
      dynamic_cast<TClonesArray*>(EvtGen.ImportParticles());
  
  vector<pythpart> hadrons;
  hadrons.clear();

  // Get the generated event information

//  float process_id = EvtGen.GetMSTI(1);
//  float parton_x1 = EvtGen.GetPARI(33);
//  float parton_x2 = EvtGen.GetPARI(34);
//  float q_squared = EvtGen.GetPARI(22);
//  float parton_pt = EvtGen.GetPARI(17);

  // Rotate the event by a phi angle
  // (pythia generates scattered partons at fixed phi)
    
  double phi_rotate = rang->Uniform()*TMath::TwoPi(); 

  Int_t numParticles = particleArray->GetEntries();
  for (Int_t i=0; i<numParticles; i++) {

      // Get particle information
      TMCParticle *particle =
         dynamic_cast<TMCParticle*>(particleArray->At(i));

      Float_t energy     = particle->GetEnergy();
      //Int_t   firstchild = particle->GetFirstChild()-1;
      Int_t   type       = particle->GetKF();
      Int_t   status     = particle->GetKS();
      //Int_t   lastchild  = particle->GetLastChild()-1;
      //Int_t   parent     = particle->GetParent()-1;
      Float_t px         = particle->GetPx();
      Float_t py         = particle->GetPy();
      Float_t pz         = particle->GetPz();

      // Event plane rotation:
       
      Float_t pyth_pt = sqrt(px*px + py*py); 
      Float_t pyth_phi = atan2(py,px);

      px = pyth_pt*cos(pyth_phi + phi_rotate);
      py = pyth_pt*sin(pyth_phi + phi_rotate);

      //cout << status << " " << type << " " << px << " " << py << " " << pz << endl; 

      if(status==1){

	//cout << type << " " << px << " " << py << " " << pz << endl; 

	pythpart temp; 

	temp.px = px; 
	temp.py = py; 
	temp.pz = pz; 
	temp.e = energy; 
	temp.type = type; 

	hadrons.push_back(temp); 
	
      }

  }

  // Write out the event

  FILE *pFile = fopen ("soft_event.txt","w+");

  //cout << hadrons.size() << endl; 
  fprintf(pFile,"%u \n", unsigned(hadrons.size()));
  for(unsigned i=0; i<hadrons.size(); i++){

    //cout << hadrons[i].type << " " << hadrons[i].px << " " << hadrons[i].py << " " << hadrons[i].pz << " " << hadrons[i].e << endl; 
    fprintf(pFile,"%d %f %f %f %f \n",hadrons[i].type,hadrons[i].px,hadrons[i].py,hadrons[i].pz,hadrons[i].e);

  }

  fclose(pFile); 



}
